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ABSTRACT 

We present time-dependent numerical hydrodynamic models of line-driven 
accretion disk winds in cataclysmic variable systems and calculate wind mass-loss 
rates and terminal velocities. The models are 2.5-dimensional, include an energy 
balance condition with radiative heating and cooling processes, and includes local 
ionization equilibrium introducing time dependence and spatial dependence on 
the line radiation force parameters. The radiation field is assumed to originate in 
an optically thick accretion disk. Wind ion populations are calculated under the 
assumption that local ionization equilibrium is determined by photoionization 
and radiative recombination, similar to a photoionized nebula. We find a steady 
wind flowing from the accretion disk. Radiative heating tends to maintain the 
temperature in the higher density wind regions near the disk surface, rather than 
cooling adiabatically. For a disk luminosity L disk = L Q , white dwarf mass M wd = 
0.6M©, and white dwarf radii Rwd = O.Oli?©, we obtain a wind mass-loss rate 
of M W in d = 4 x lO~ 12 M 0 yr _1 , and a terminal velocity of ~ 3000 km s J . These 
results confirm the general velocity and density structures found in our earlier 
constant ionization equilibrium adiabatic CV wind models. Further we establish 
here 2.5D numerical models that can be extended to QSO/AGN winds where, the 
local ionization equilibrium will play a crucial role in the overall dynamics. 
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Subject headings: accretion, accretion disks — hydrodynamics — novae, cata- 
clysmic variables — stars: mass-loss 


1. Introduction 

Ultraviolet observations of Cataclysmic Variables (CVs) show that virtually all non- 
magnetic CVs with high mass accretion rates (>4x lO -lo M 0 yr -1 ) exhibit P Cygni profiles 
in the resonance lines (e.g. Cordova k Mason 1982; Greenstein k Oke 1982; Guinan k 
Sion 1982; Klare et al. 1982; Vitello k Shlosman 1993; Prinja k Rosen 1995; Friedjung, 
Selvelli, k Cassatella 1997; Gansicke, Beuermann, k Thomas 1997; Knigge et al. 1997; 
Prinja et al. 2000). These profiles consist of blue-shifted absorption and redshifted emis- 
sion, with the emission most prominent in high inclination systems (e.g. Hutchings 1980; 
Krautter et al. 1981; Holm, Panek, k Shiffer III 1982; Cordova k Mason 1985; Mason 
et al. 1995; Cheng et al. 1997), and the absorption appearing at low and intermediate 
inclinations. The most straightforward explanation for these observations is in terms of a 
wind originating from the accretion disk. The wind is driven by radiation pressure from the 
disk UV continuum absorbed by line transitions, analogous to the radiation driven winds in 
early type stars. Support for this scenario came from Drew (1987), who calculated theoret- 
ical C IV 1549 A line-profiles from one-dimensional kinematic CV wind models, for a wide 
variety of wind and disk parameters. Drew (1987) found that the dependence on the form 
of the UV resonance lines with inclination angle could be accounted for by a bipolar wind 
emerging from an accretion disk. 

One-dimensional dynamical models for line-driven winds (LDW) in CVs were developed 
independently by Vitello k Shlosman (1988) and Kallman (1988). Further efforts came 
through two-dimensional kinematic modeling (Shlosman k Vitello 1993; Knigge, Woods, k 
Drew 1995), which succeeded in showing consistency between the assumed polar geometry of 
a disk wind and observed profiles; but such models were unable to self-consistently calculate 
the wind dynamics including the effects of rotation and the anisotropic disk radiation field. 

Two-dimensional disk wind dynamical models were developed by Icke (1980), but these 
did not take into account the radiation pressure due to line absorption. For a typical white 
dwarf mass of O.6M 0 (Leibert 1980; Silvestri et al. 2001) the disk luminosity required 
to produce such a wind without line radiation pressure (assuming radiation pressure due 
to continuum scattering only) would require a luminosity of ~ 1O 5 L 0 (Icke 1980), several 
magnitudes above the observed luminosity of CVs which vary between 0.01 L 0 and 1OL 0 
(Patterson 1984). Icke (1980) did find that, for sufficiently high disk luminosities, a biconical 
disk wind would form. Icke (1981) suggested that biconical winds were a general property 
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of accretion disk winds independent of the wind driving mechanism. 


Two-dimensional dynamical models of isothermal line-driven disk winds were presented 
by ns Pereyra, Kallman, k Blondin (1997, hereafter Paper 1). Results from Paper 1 showed, 
in analogy with line-driven winds from early type stars, that terminal velocities are approxi- 
mately independent of the luminosity of the disk, although increments in luminosity produce 
increments in mass-loss rate, and that rotational forces are important in the study of winds 
from accretion disks. They cause the velocity streamlines to collide, which reduce the speed 
and increase the density of the wind producing an enhanced density region The highes 
absorption occurs in the enhanced density region where density is increased relative to a 
spherically diverging wind with the same mass loss rate and the velocity is roughly half the 
terminal velocity. This density increase is necessary in order to produce at least margina y 

optically thick lines. 

We developed (Pereyra 1997; Pereyra, Kallman, k Blondin 2000, hereafter Paper 2) an 
adiabatic wind model (rather than isothermal), driven by a standard accretion disk (Sha ura 
& Sunyaev 1973; Lynden-Bell k Pringle 1974) rather than an isothermal disk. The hydro y- 
namic models developed of line-driven accretion disk winds (LDADWs) m CVs me u e e 
radial structure of an optically thick accretion disk with the corresponding radiation fields 
and surface temperature distributions. The corresponding energy conservation equation was 
also implemented, including the adiabatic heating and cooling effects due to compression 
and expansion. From the computational models, assuming single scattering and constant 
ionization, we calculated theoretical line profiles for the C IV 1550 A line, and found that 
the line profiles obtained were consistent with observations in their general form and strong 
dependence with inclination angle. 


Proga, Stone, k Drew (1998) also developed two-dimensional models for these systems, 
obtaining similar results to the models presented here. They assumed an isothermal wind 
and constant values for the line radiation parameters throughout the wind. They found 
similar wind velocities and wind mass loss rates as we do here, although they assumed dis 
luminosities about an order of magnitude greater. A difference we find with the results of 
the models presented by Proga, Stone, & Drew (1998) is that they obtain unsteady flows 
characterized by large amplitude fluctuations in velocity and density, while for our previous 
models (Paper 1; Pereyra 1997; Paper 2) and in the models presented here we find a steady 
flow. We note that these large fluctuations (rather than the steady ow we n ) cou 
account for the difference in luminosity required by the Proga models, with respect to our 
models, to obtain similar wind mass loss rates. One possibility suggested by Proga, Stone, 
k Drew (1998) was that the spatial resolution of the models we presented in Paper 1 was 
too low to adequately resolve the wind. In paper 2 we increased significantly the spatial 
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resolution of our models and obtained similar steady results when compared to our earlier 
spatial resolutions used in Paper 1; thus showing that differences in spatial resolution of 
t e models was not generating the difference between our steady disk wind models and the 
intrinsically unsteady” disk winds reported by Proga, Stone, k Drew (1998). Another 
possibility is that the difference in the unsteady flows of the Proga models with respect to 
t e steady flows we find is due to the boundary condition treatment in the Proga models 

which appears to lead to instabilities at the base of the wind and in turn to the strong 
fluctuations reported. 6 

It is well known that for early type stars, when a Sobolev treatment is used for the 
hne radiation pressure (as we have done in our models and as has been applied by Proga 
Stone & Drew (1998) for accretion disk winds in CVs), a steady wind solution is obtained 
(e.g. wocki k Puls 1999). The treatment of boundary conditions in numerical models for 
me driven winds in early type stars was discussed by Owocki, Cranmer, k Blondin (1994) 
and Cranmer k Owocki (1996). They find that the treatment of the lower (base of wind) 
boundary condition is “somewhat problematic”. Further they find that if they used an 
arbitrarily high density, similar to the boundary conditions used by Proga, Stone k Drew 
(1998), they would find oscillations at the base of the wind in both density and velocity. 
In both papers, in addition to other details in the boundary condition treatment (Owocki. 
Cranmer, k Blondin 1994; Cranmer k Owocki 1996), they lowered the density at the inner 
oun arj until they found steady conditions at the base and were simultaneously able to 
resolve acceleration in the subsonic region. Therefore a conclusion in the area of numerical 
modeling that we derive is that by varying the treatment of boundary conditions at the base 
of a line driven wind, one can obtain strong fluctuations in density and velocity which are 

not necessarily physically intrinsic but are rather a consequence of the numerical details of 
the model. 

Feldmeier k Shlosman (1999) and Feldmeier, Shlosman, k Vitello (1999) developed 
two-dimensional stationary models for CV disk winds where they solve the dynamics under 
the assumption that the wind streamlines lie on straight cones. The values of wind mass 
loss rates obtained through the models presented by Feldmeier, Shlosman, k Vitello (1999) 
for an isothermal disk (equation [11] of their paper) are consistent with those found in 
Paper 1. In addition Feldmeier, Shlosman, k Vitello (1999) find that in an isothermal 
disk (equations [10]-[11] of their paper) the wind mass loss rate scales as M wind oc L l J° k , 
where L disk is the disk luminosity and a is line force multiplier parameter. This is the 
same scaling law reported in Paper 1. Further, the biconical geometry and wind tilt angle 
in the models of Feldmeier, Shlosman, k Vitello (1999) are roughly consistent with the 
velocity structures we find through full 2.5D hydrodynamic modeling of our earlier models 
(Paper 1; Pereyra 1997; Paper 2) and in the models presented here. We further note that 
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Feldmeier et al. suggest that detailed non-LTE calculations of the line-force multiplier may 
be important in the modeling of CV disk winds. This work is an additional step in that 

direction. 

In spite of these encouraging results, a more realistic model which would include a 
detailed calculation of the UV spectrum from dynamically calculated density and velocity 
fields is desirable in order to make specific comparisons with observations. For instance, 
a maximum flux greater than 1.7 times the continuum has been reported for the emission 
component of P Cygni line profiles of C IV 1550A from some cataclysmic variables (e.g. 
Cordova & Mason 1982), while from our previous models (Pereyra 1997, Paper ) we n 
that the emission maximum of the C IV lines stays below 1.3 times the continuum flux for 
those angles where the absorption component is observed. This discrepancy may be due to 
some of the simplifications in our previous models or it may indicate a C IV line-emitting 
region other than the disk wind itself. Furthermore, in recent observational results, Froning 
et al (2001) find relatively narrow absorption lines (< 700 km s ) in the FUV spectra o 
the cataclysmic variable U Gem. As Froning et al. (2001) have indicated, such lines cannot 
arise from the disk photosphere. In their work they suggest the possibility of an outer disk 
chromosphere to account for these lines, but it is also possible that some of these narrow 
FUV lines could be accounted for by an accretion disk wind. 

An additional motivation for this work is that the models developed here can also be 
extended to the study of LDADWs in Quasars (QSOs), where the local ionization equ.hbnum 
will play an important in the overall dynamics (Murray et al. 1995; HUUer et al. 2002), 
Evidence for winds in QSOs is found in the broad absorpt.on lines (BALs) observed in 
approximately 10% of the QSOs (e.g. Weymann 1997). If the wind in QSOs is generated 
from an accretion disk, then the presence of BALs in only a fraction of QSOs can be accounted 
for as a viewing-angle effect (Tnmshek 1984); similar to the case of CVs where the accretion 
disk wind produces P Cygni profiles when observed “face-on” and emission lines ■ (without 
absorption components) when seen “edge-on" (Pereyra, Kallman, k Blond, n 2000). The 
x-ray luminosity in QSOs are generally comparable to the UV/optical luminosities (eg. 
Tananbaum et al. 1979; Grindlay et al. 1980; Mushotzsky, Done k Pounds 1993, 
Laor 1997; George et al. 2000). Drew k Boksenberg (1984) found that such high x- 
ray luminosities would ionize the wind to a point where the populations of ions responsible 
for the BAL observed in QSOs would be too low to produce observable absorption lines. 
Murray et al. (1995) found that, with an appropriate x-ray shielding mechanism LDA s 
can produce wind densities and velocities consistent with BALs observed in QSOs^ Strong 
changes in the ionization balance can produce observable effects in line formation in QSOs; as 
well as significant changes in the line radiation force parameters which will play an important 
role in the overall wind dynamics (Hillier et al. 2002). 
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In our previous models (Pereyra 1997; Paper 2) we had assumed constant ionization 
equilibrium throughout the wind and single scattering in the line profile calculations. In this 
work we extend our previous models to include radiative heating and cooling and ionization 
balance in the wind, and a calculation of the line radiation force parameters at each point 
in space at each time step. Our results are similar to those of Paper 2, and therefore act to 
confirm the general results found in our previous simpler models. They are a further step 
towards a model capable of self-consistently producing theoretical spectrum which may be 
compared in detail with observations; they can also be extended to the study of LDADWs 
in QSOs/AGN, where the local ionization equilibrium will play a crucial role in the overall 
dynamics (Murray et al. 1995; Hillier et al. 2002). 

In §2 we discuss the radial structure of the accretion disk and the radiation field as 
implemented in our models. In §3 we present the ionization balance calculations applied in 
this work. We derive the expression used for the treatment of the line radiation pressure in 
§4. The treatment of radiative heating and cooling in this work is discussed in §5. In §6 we 
present and discuss the hydrodynamic calculations. The results of the models are presented 
and discussed in §7. We present a summary and the conclusions of this work in §8. 


2. Accretion Disk 

A key ingredient in the models presented in this work is the radiation field emerging 
from the accretion disk. The boundary layer where the disk intercepts the white dwarf and 
the white dwarf photosphere may also contribute to the radiation field. The white dwarf by 
itself (without considering the effects of the accreting mass) will have a luminosity of the 
order of O.OILq (Leibert 1980), while the accretion disk will have typical luminosities of the 
order of L e (Warner 1987). As we indicated in Paper 2, there is still uncertainty concerning 
the existence and spectrum of the boundary layer (Vrtilek et al. 1994; Mauche, Raymond, 
&; Mattei 1995, Long et al. 1996; Szkody et al. 1999; Mauche k Raymond 2000; Froning 
et al. 2001). In the models presented here we have chosen not to explore the boundary 
lajer issue in depth, but instead to consider the non-boundary layer scenario. We assume 
m this work that the radiation field is generated by the accretion disk alone, neglecting 
any boundary layer or white dwarf radiation. We assume here the standard accretion disk 
(Shakura k Sunyaev 1973; Lynden-Bell k Pringle 1974). In paper 2 we discussed some of 
the observational evidence for the existence of an accretion disk in CVs. We note that in 

this work we study the winds originating in nonmagnetic CVs in which an accretion disk is 
formed. 

In the standard accretion disk model, the disk is assumed to be in a steady state. Shear 
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stresses transport angular momentum outwards as the material of the gas spirals inwards. 
Conservation of angular momentum leads to the following expression: 

W2nr 2 - M a ur 2 = C (constant) (!) 

where r is the radius, M a is the mass accretion rate, u is the corresponding angular velocity, 
and W is defined by: 

W= [ Wrd> dz ( 2 ) 


-20 


where z 0 is the half thickness of the disk and w r<p is shear stress between adjacent layers. 
With the additional assumption that the shear stresses are negligible in the inner disk ra ms: 

Ma 


W=1 p-J UJ r 2 -uJor 2 ) 
2'kt 1 


( 3 ) 


where u, 0 is the angular velocity at the inner disk radius r„. In CVs the inner disk radius 
corresponds to the white dwarf radius R wd , i.e. the accretion disk extends from the inner 
Lagrangian point of the binary system down to the white dwarf surface. 

Taking into account the work done by shear stresses, and assuming that, as the mass 
accretes inwards, the gravitational energy lost is emitted locally: 


3=737 


1 d 


47 rr dr 


. fco 2 r 2 GM\ 
Ma V 2 r ) " 
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( 4 ) 


where Q is the radiated energy per area of the disk surface. Assuming that the disk material 
follows Kepler ian orbits, i.e. 

u = (5) 

where G is the gravitational constant and M wd is the white dwarf mass; from equations (3)- 
(4) one finds: 


Q = 


3 M aC cr G M w d 
87rr 3 


[Ht) 


1/2 


( 6 ) 


The function Q{r) (equation [6]) was originally derived by Shakura k Sunyaev 197 
for binary systems with an accretion disk about a black hole. Popham k Narayan (199 ) 
developed detailed models of boundary layers and accretion disks in cataclysmic variables 
and found that their models supported the same function Q(r) for the accretion disk in CVs. 

Further, assuming that the disk is emitting locally as a blackbody, the radial temperature 
distribution of the disk will be given by. 

x 1/4 
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3 Maccr G M w d 


87rr 3 cr s 


( 7 ) 
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where <r s is the Stefan-Boltzmann constant. Equation (6) is implemented in our model in 
the calculation of radiation flux throughout the wind, and consequently in the calculation of 
the continuum and line radiation pressure force. We assume local blackbody radiation. In 
our model we also assume that the temperature at the base of the wind is equal to the disk 
surface temperature for each radius. We wish to note here that the blackbody assumption of 
equation 7 will not strictly hold. Radiative diffusion of the energy from the inner parts of the 
disk to the disk photosphere will produce deviations from the blackbody assumption, limb 
darkening effects will produce deviations from the projected area in the angular dependence 
of the emitted flux, line scattering and electron scattering will produce deviations from the 
blackbody emission, the relatively high velocity gradients in the rotational shears as well as 
the wind velocity gradients will affect the line scattering; thus a detailed calculation of the 
spectrum emitted from the accretion disk would requires a 2.5D radiative transfer calculation 
of the disk atmosphere coupled with the corresponding velocity gradients. 

Early attempts to calculate deviations of the actual disk spectrum from the local- 
blackbody assumption, applied locally spectra from stellar atmosphere codes (Herter et al. 
1979; Kiplinger 1979, 1980; Mayo, Wickramasinghe, k Whelan 1980; Kolykhalov k Sun- 
yaev 1984). More recently, local plane parallel radiative transfer models implementing disk 
vertical structure have been developed (Diaz, Wade, k Hubeny 1996; Hubeny k Hubeny 
1997, 1998; Wade k Hubeny 1998; Hubeny et al. 2000). However, to our knowledge a self- 
consistent 2.5D radiative transfer model coupled with velocity shears of an accretion disk 
atmosphere has not yet been developed. Also, there still remains uncertainties regarding 
the variation of viscosity in the disk which in turn effects the vertical structure. Thus, due 
to its simplicity, the assumption of local blackbody spectrum has often been used to model 
accretion disk and their spectra (e.g. Shlosman k Vitello 1993; Knigge, Woods, k Drew 
1995; Proga, Stone, k Drew 1998; Pereyra, Kallman, k Blondin 2000). In future work 
we shall consider the radiative transfer effects that deviate the disk spectrum from an ideal 
local-blackbody, and attempt to couple them with hydrodynamic calculations. 

It is generally accepted that the mechanism for producing the shear stress in accretion 
disks is magnetic turbulence generated by the Balbus-Hawley magnetorotational instability 
(MRI) (Balbus k Hawley 1991). Three-dimensional magnetohydrodynamic simulations (e.g. 
Armitage 1998; Machida, Hayashi, k Matsumoto 2000; Hawley, Balbus, k Stone 2001; 
Hawley k Krolick 2001) suggest that the standard assumptions (e.g. Keplerian orbits; 
zero stress at inner radii) may not hold and that therefore corrections to equation (6) may 
arise. We believe that the standard Shakura-Sunyaev disk includes much of the physics 
involved in accretion disks, such as: gravity from a large central mass determining the 
aveiage velocity fields in the disk, angular momentum transport, conversion of gravitational 
energy to radiative emission, and monotonic increase of temperature with decreasing radius 
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(except of radius close to the inner disk radii) and local blackbody emission which accounts 
for the relatively “flat” continuum spectrum observed in systems where accretion disk are 
inferred. We also believe that it is important to explore “simple” models where the number 
of free parameters is small. More realistic disk models and radiation fields will be introduced 

at a later date. 


3. Ionization Equilibrium and Line Opacity 

The main wind driving mechanism in our model is the line radiation force generated 
as the disk emission encounters ionized gas above the disk. The wind material is thus 
driven from the disk photosphere achieving velocities of the order of a few thousand of 

kms -1 The line radiation force will depend on the available line opacity distribution o t e 

wind material, which in turn is a direct function of the ion populations. Therefore the ion 
populations become an important ingredient in our models. In our previous work ( aper , 
Pereyra 1997, Paper 2), in an approach similar to that presented by Castor, Abbott, & Klein 
(1975) for OB star winds, we had assumed constant ionization equilibrium throughout the 
wind. Further we had assumed that the ion populations were similar to those of earl} 
type stars. These two assumptions allowed us to apply typical OB star line radiation force 
parameters (Abbott 1982) throughout the wind, significantly simplifying the numerica 
treatment. In this work we include the local ion population calculations simultaneously with 
the hydrodynamics as we describe below. 

The ion populations of the wind material is initially at local thermodynamic equilibrium 
(LTE) since the wind starts at the disk photosphere. As the wind material rises, the ion 
populations are no longer in LTE. The material is irradiated by the disk, which is a non-Plank 
source (equations [6]-[7]). Due to the relatively low wind density, the ion populations above 
the disk are largely determined by photoionization and radiative recombination processes, 
which in turn are dependent on the incoming disk radiation. It is difficult to model t e 
collisional and radiative processes which determine the ion populations, varying in space and 
time, while simultaneously solving the hydrodynamics of the system. In this work we assume 
that the ion populations in the disk wind are determined by photoiomzation.and radiative 
recombination processes, similar to a photoionized nebula irradiated by an accretion is • o 
illustrate the validity of this approximation we shall estimate the ratio between the ra lative 
photoionization rate Rtj and the collisional ionization rate C iy Under the condition la 
hv » kT R and hu 0 » kT e , where h is Plank’s constant, hv Q is the minimum energy require 
to ionize an atom from state i to state j, k is Boltzmann’s constant, T R is the temperature 
of the ionizing radiation, and T e the electron temperature, the ratio can be approximated by 
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(Mihalas 1978): 


Rjj ~ 4(2tt 3 A;)U^ ( WTr 


c, 


y 


3m e 2 e 2 c 3 


ji e T) 


exp 


M-L- Ml 


(8) 


where m e is the electron mass, e the electron charge, c the speed of light, W the dilution 
factor of the radiation from its source, and n e is the electron density. Taking for the CV 
case. n e ~ 10 (this density is a typical value for the CV disk wind models we presented in 
Paper 2), T R ~ T e ~ 30, 000 K , kT ~ 1/2, and hv 0 ~ 13.6 ev, from equation (8) we find: 


^~1.4x 10 8 . 


( 9 ) 


As we have indicated above, in order to calculate the force due to line scattering, one 
in principle must consider the contribution to all lines which in turn will depend on the ion 
populations at a particular point in space in a particular time. In the calculation of ion 
population and line opacity we follow a similar scheme as that applied by Abbott (1982) for 
OB stars. In this work we assume that photoionization and radiative recombination processes 
determine the ion populations, therefore the ionization balance in the wind is determined 
through Saha s equation for a gas photoionized by an external source: 


Aj-f l Tl e 


2IP^±i 

Uj 


2n m e kT e 
h 2 


) 


3/2 




(10) 


where N 3 and N j+l are the atom densities at stage j and at the following ionization stage 
j + 1, Uj and U J+l are the partition functions for the corresponding ionization stages, and 
Xij is the ionization potential for the atom at stage j. Equation (10) has been applied to 
several astrophysical systems, e.g. stellar winds (Abbott & Lucy 1985; Schmutz et al. 1990- 
Lucy & Abbott 1993). 


In the work presented here, we have assumed that the radiation is primarily coming from 
the accretion disk which presents a temperature dependent on radius (equation [7]). Equa- 
tion (10) assumes that the photoionization rates are much greater than collisional ionization 
rates, thus: 

Njlj « Nj + iRj + i ( 11 ) 

where Ij is photoionization rate for an atom at ionization stage j , and R j+l is the radiative 
recombination rate for the atom at ionization state j + 1. Therefore, from equation (10): 


*i = R 


■ 7 + 1 


2lVU J+l (2 nm e kT e \ 3/2 T R 
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Xjj 


( 12 ) 
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In the accretion disk, the radiation is emitted at different temperatures (equation [7]), thus 
the total photoionization rate I Tj due to the disk radiation must equal the sum of all the 
contributions of photoionization rates due to the radiation emitted at each section of t e 

disk surface: 


hj 


— i?j+i 


2 Uj + 1 

n e Uj 



where W{r,<j)) is the dilution factor per source area, given by: 


Xi± 

kT R (r) 


r dcj) dr 


(13) 


W{r, </>) = 


1 z- (fp-r) 
47T I fp - r\ 



(14) 


ie - : 1 Z Clbl 

w(r ’ V = 4i [IP - 2Rr cos + r> + W 2 

where r p 3 (R, 0, z) is the position where the dilution factor is evaluated and r = (r cos 0. r sin <fr, 0) 
is the radiation source position on the disk. 

Therefore, under the assumption that the photoionization and radiative recombination 
processes dominate over the collisional processes, we have: 


(“) 3 / 2 sL^™ exp (-^)) r#iir (16) 

which is the equation that we implement in this work for the calculation of ionization balance 
for the line radiation force. 


Once the ion populations at different stages is calculated for each atom through the 
resolution of equation (16), the population for each energy level for a given ion stage of a 
given atom is obtained through Boltzmann’s equations: 


Nil = 

N 0J 


9ij 

9oj 


exp 



(17) 


where N„ and N oj are the densities of a given atom at a given ionization stage j at energy 
level i and at the lowest energy level respectively, ft, and ()„, are the degeneracies of energy 
level . and at the lowest energy level respectively, and X tj is the difference in energy between 
level i and the lowest energy level. We wish to note here that the LTE excited population 
assumption in equation (17), while allowing the local ionization equilibrium calculations to 
stay numerically manageable, will tend to be more accurate for the lower energy -on states 
which are the most crucial for the line radiation force. For early type stars Abbott (1982) 
found that roughly 90% of the line radiation force was due to lines arising from groun or 

metastable states. 
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After the ion populations at different stages at different energy levels for a given atom 

have been obtained, the line opacity of a given line for a given atom is calculated by the 
expression: 






1 

A v D 


(18) 


where g h and g UJ are the degeneracies of the lower energy level l and the upper energy 
level u of ionization stage j for a given line of a given atom, f luj is the absorption oscillator 
strength corresponding to the transition between lower energy level l and upper energy level 
U of i°nization stage j for a given line of a given atom, N tj and N uj are the ion densities of 
ionization stage j, lower energy level l, and upper energy level u respectively, p is the total 

mass density of the gas, and /\u D is the Doppler width of the given line. The Doppler width 
of the line is given by: 


Ai' 



(19) 


where v th is the thermal velocity and the line frequency u L is given by: 




X U] - X tj 
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( 20 ) 


4. Radiation Force 


Once the line opacity k l has been calculated through the above expressions for each 
atomic line at each ionization stage, the line radiation force is given by: 
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( 21 ) 


where I{h) is the radiation intensity, and r L (n) is the line optical depth. 

The line optical depth r L (n) for a given line at a given direction n, under the Sobolev 
approximation, is given by (Rybicki k Hummer 1978): 


r L (h) 


P^LVth 
n ■ V(n ■ wj| 


( 22 ) 


As in Paper 1 and Paper 2, we assume in this work that the velocity gradient is primarily 
in the “z” direction and approximate the line optical depth by: 

n(h) * - gat . 

| ft ■ z dv z /dz\ 


( 23 ) 
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Integrating equation (21) in optical depth. 


lines 


1 - exp (-r L (n)) 

n(n) 


1(h) hd£l 


(24) 
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where cr is Thomson cross section per unit mass. 

In their work on radiation driven winds in OB stars, Castor, Abbott, & Klein (1975) 
and Abbott (1982), found that the sum over lines in equation (25) could be approximate 

by the following expression: 


1(h) Aud 
f n °° I(h)dv 

lines ■'O 
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ex P ( - f i^ T r(") ) 


^r t (n) 


fe'(n) (— 

\KL 


n(n) 


- a'(n) 


(26) 


where the k'(h) and a'(h) are constants for a given direction n. Defining t: 


t = — r L (h ) 

KL 


we have: 


E 


I (h)Avo 



k'(h)t' 


-o'(n) 


(27) 


(28) 


lines l T 1 ^ , 

From equation (28) we obtain the fc'(n) and the a'(h) in a form similar to the Abbott 
(1982) calculations for OB stars. We calculate the left hand side of equation (28) varying 
the t parameter from KT 1 to KT 7 , using the line opacity <c t calculated in the form describe 
in the previous section for each line, and using the disk emission described by equations (6)- 
(7) to calculate the intensities. From these values, the parameters k (n) and a (n are 
obtained through linear regression by varying the t parameter over two order o magni u es 
at a time within the fore mentioned interval, and taking the highest values foun or ea 
parameter. We note here that the selection of the highest values within intervals of two orders 
magnitudes is somewhat arbitrary. To our knowledge this is the first attempt to me u e e 
spatial and time dependence of the line radiation force parameters coupled consistently with 
the hydrodynamics of accretion disk winds, and we believe that the process includes much 
of the physics involved in line driving from a distributed radiation source: line-scattering, 
photon momentum generating force on gas, dilution factor, line optical depth in a moving 
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media, line-opacity distribution, energy distribution in spectrum, and ion populations This 
approach allows us to explore the effect of variable ion populations on the disk wind rather 
that assuming a constant ionization equilibrium as we did in Paper 2. In future work we 
shall develop a more accurate treatment for the line radiation force parameters, including a 
systematic determination of the appropriate intervals in t in which to derive the line force 
parameters, as well as a study of the directional dependence of these parameters, rather than 
assume them isotropic as we do in this work (see equation [30]). 

From equation (25) and equation (28), we have: 




/(n)Gfc/ ^ 

c 


( h)t 


nd£l . 


(29) 


s mentioned above, in this work we use values of the force multiplier parameters k and 
a which vary in time and space but are independent of direction. As we see below this will 
a ow us to implement the numerical integrations on intensity 1(h) at each spatial point of the 
computational grid before implementing the hydrodynamic calculations rather than during 
the implementation of the hydrodynamic calculations, thus maintaining the computations 


The values of k and a used in this work are calculated through the expression 


/ 

f 0 c ° I(h)dvk'(h)t~ a 'W 

hdQ, 

1/ [Jr I{n)dv] hdCl 



e calculate the right hand side of equation (30) varying the t parameter from 10' 1 to 
, using the values of k'(h ) and a'(h ) calculated through equation (28) and using the 
disk emission described by equations (6)-(7) to calculate the intensities. From these values, 
the parameters k and a are obtained through linear regression. We note here that the angle 
averaging is done in an ad hoc manner in order to reduce computational complexity. As we 
indicated above, we believe that this approach includes much of the physics of line driving 
within an accretion disk wind. In future work, we shall develop a more systematic treatment 
of the line radiation force parameters that will include a study of the directional dependence 
of these parameters, rather than assume averaged angle values as we do here. 

With these values of k and a, we approximate the line radiation force by: 



zT [^a kt 

c 


hdQ, . 


( 31 ) 
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Therefore, from equation (27): 

- f \ a r 0 °°/(n)^ fa 
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Substituting equation (23): 
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Thus: 


= 1 




dv z 
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In • z\ a / I{n)du 

dz 

) j 
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71 dQ . 


(34) 


We wish to note that the k and a parameters in the above equation depend on the line 
opacity distribution (equations [26]-[30l), which in turn depends on density and temperature 
(equations [15]-[18]). Since the density and temperature will vary in space and time e 
k and a parameters will also depend on space and time. Also, we wish to note here that 
in this work we assume azimuthal symmetry and thus have only two independent spatial 

coordinates “r” and “ 2 ”. 

In equation (34) we substitute the a parameter inside the integral by a fiducial value of 
0 7 which is the value we assumed in our previous models presented in Paper 1 and Paper 2. 
Varying the number density from KPcnr 3 to 10 10 cm- 3 at spatial points throughout our 
computational grid, and applying the values of disk luminosity L disk = L Q and white dwarf 
radii of R, = 0 OIK© which we use throughout this work, by the way of equation (30) 
« And values for the o parameter which vary from 0.55 to 0.85. We do not expect that 
the application of the fiducial value for a inside the integral of equation (34) will produce a 
significant difference in the density and velocity wind structures, and furt er it wi a ow e 
evaluation of the integral at all spatial grid points before solving the hydrodynamic equations 
described in §6, maintaining the computations manageable. Thus, 


frad-^Hr,z,t) 


( 1 
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\pcrvth 
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(35) 


We wish to also note that although equation (35) takes into account the line radiation 
pressure in the V’ direction, we have included only the derivative m the “z” direction of he 
‘V component. As we found in Paper 2, the accretion disk winds in CVs tend to How _,n it 
direction perpendicular to the disk, particularly in the region near the surface °f ‘he disk, 
justifying the assumption that the derivatives of the ‘V’ component are neg igi e. 
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work we also assume negligible the radial derivatives of the ‘V’ component and the vertical 
and radial derivatives of the %” component in the calculation of the line radiation force. 
We note here that the disk wind streamlines are “helical” in nature; the wind flow starts 
out at the base of the wind with a non-zero component due to disk rotation. The “tu” 
component of the velocity fields is calculated throughout our hydrodynamic calculations 
under the assumption of azimuthal symmetry (see Appendix). We note that even under 
azimuthal symmetry, line radiation force may have non-zero components along the the “0” 
direction; e.g., Gayley & Owocki (2000) for the case of OB stars, assuming azimuthal 
symmetry, showed that the 4> component of the line radiation force, although presenting 
modest modifications in the overall wind dynamics, could lead to observable effects. Even 
though, on one hand, the inclusion of the additional terms in the velocity gradient for the 
line radiation pressure in the accretion disk wind may possibly produce interesting results, 
on the other we do not believe that it will significantly affect the overall results currently 
being implemented. In future work we shall include the effects of the derivatives along both 
the “ 2 ” and “r” directions for the different velocity components in the line radiation force, 

rather than assume “ dvjdz ” as the dominant term as we have done in this work and in 
Paper 1 and Paper 2. 

Defining S(r, z ): 


S(r, z) 


\fi • i| 0,7 I(n)du 


hdQ 


(36) 


Taking the source of radiation to be the accretion disk, we have: 

O n 27T 


,(r,z)= r f 
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Jr Q JO 


7T f(r 2 + r' 2 + 2 2 - 2rr , cos0) 1 /2]4.7 r d( ^ dr 


(37) 


(38) 


where S z and S r are the corresponding components of the vector S defined in equation (36) 
(the 5’ 0 component is zero due to the axial symmetry of the accretion disk) and Q(r') is the 
radiation emission per area of the disk calculated through equation (6). 


Thus we have 


frad — k(r, Z, t) ( 

c V pcrv th 


dv. 


dz 


a(r,z,t) 


[Sz(r, z)z + S r (r,z)f] 


(39) 


where S z and S r are given by equations (37) and (38). 
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We also account for the fact that, as discussed by Castor, Abbott, & Klein (197o) for 
OB star winds, if the velocity gradient increases to sufficiently high values, or the density of 
the wind decreases to sufficiently low values, the contribution for the radiation pressure due 
to lines arrives at a maximum value. Similar to the case of OB stars (Abbott 1982) we found 
that the maximum value of the line radiation pressure was obtained when (1/ pav th ) dv z /dz „ 
10 7 In this work we have assumed that the maximum possible value of the line radiation 
pressure is obtained when (1 /pav th ) dvjdz = 10 7 (rather than 10 8 as we did in Paper 2). 

Therefore the expression for the total line radiation force per mass for disk winds which 
we adopt is 


frad = ^Hr, z,t) (max 


1 

J pOVth 


dvz 

dz 


,10 7 


a{r,z,t) 


[5 z (r,z)z + 5 r (r,z)r] . (40) 


In order to calculate the k(r,z,t) and a(r,z,t) functions self consistently throughout 
the hydrodynamic calculations described below and simultaneously keep the computation 
manageable the numerical code was written such that the spatial distribution of the k an 
a are calculated for values of number density varying from 10 6 cm" to 10 cm at each 
computational spatial grid point. A fiducial temperature is taken for these calculations, since 
as we find in the models presented here, in inner disk region where the line radiation force 
is most crucial, the temperature does not change significantly due to radiative heating an 
cooling During the hydrodynamic calculations, the actual values of k and a are obtaine 
through geometric interpolation from the above values and the current value of density at 

each computational grid point. 

For the continuum radiation pressure the total radiation flux is calculated throughout 


the disk wind: 



hdQ. . 


Taking the source of radiation to be the accretion disk, we have 


(41) 
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r'd<t> dr' 
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and 


F r (r , z) 



z(r — r' cos (j>) 

[(r 2 + r' 2 + z 2 — 2rr'cos<£) 1/2 ] 4 


r'dcj) dr' 


(43) 


where F z and F T are the corresponding components of the radiation flux vector F originating 
from the disk (the component is zero due to the axial symmetry of the accretion disk) 
and Q{r') is the radiation emission per area of the disk calculated through equation ( ). 
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5. Radiative Heating and Cooling 

Radiative heating and cooling are computed in our model assuming an optically thin 
disk spectrum. A local balance between heating and cooling is assumed, subject to the con- 
straint of particle number conservation for each element. The ionization processes included 
m the calculation are photoionization, charge transfer ionization, and collision^ ionization. 
Recombination processes include radiative and dielectronic recombination and charge trans- 
fer. The heating terms includes photoionization heating, Compton heating, charge transfer, 
and colhsional de-excitation. Cooling terms include radiative and dielectronic recombination, 
bremsstrahlung, colhsional ionization, collision excitation of bound levels, and (endothermic) 
charge transfer. Details of the radiative heating and cooling calculations are described in 
detail by Kallman k McCray (1982). 

Disk models have succeeded in accounting for the observed continuum spectrum in 
many CVs (e.g. Herter et al. 1979; Kiplinger 1979, 1980; Bath, Pringle, k Whelan 1980- 
Mayo, Wickramasinghe, k Whelan 1980; Krautter et al. 1981; Hassall 1985; Wade 1988- 
la Dous 1989; Diaz, Wade, k Hubeny 1996; Wade k Hubeny 1998). The form of the 
spectrum applied in this work for the radiative heating and cooling calculations is obtained by 
integrating contributions over the disk surface calculated from equation (6), assuming local 
blackbody emission, and viewing the disk at infinity. In Figure 1 we show the disk spectrum 
used for the radiative heating and cooling calculations. We note that an assumption made 
m the radiative cooling and heating calculations is that the form of the spectrum is position 
independent, and thus the rates calculated here will be more accurate higher above the disk. 
In the future we will introduce more accurate heating and cooling rates into our models. 

In Figure 2 and Figure 3 the radiative heating and cooling are shown respectively as a 
function of temperature for different values of the ionization parameter f defined as: 

c _ 47^1 

C = — ( 44 ) 
where F is the radiation flux and n the number density. 

From Paper 2 we found that the wind temperature of CV winds range in values of up 
to ~ 40, 000 K. Thus in our work here we have approximated the heating and cooling rates 
the results presented in Figure 2 and Figure 3 by analytical expressions that agree with the 
detailed calculation within 20% for the heating rates and within 30% for the cooling rates 
for temperatures of up to 50, 000 K. 

For the radiative heating rate we implemented the following expression (in units of 



-19- 


ergcm 3 s *): 


r = 1.6 x 10 -21 T~ 1/2 


O-tot) 


(45) 


.3 x 10 5 

where T is the absolute temperature, valid for temperatures of up to 50, 000 K. 

For the radiative cooling rate we implemented the following expression (in units of 

ergCm3S " ): A ,1.9x10-7** W 

where T is the absolute temperature, valid for temperatures of up to 50, 000 K. 


6. Hydrodynamic Model 

Our hydrodynamic model uses the Piece-Wise Parabolic Method (PPM) numerical 
scheme (Colella & Woodward 1984). We use the 2.5-dimensional hydrodynamic equations 
of the wind in cylindrical coordinates; i.e., the three-dimensional hydrodynamic equations 
are reduced to two-dimensional equations by assuming that the derivative of any physical 
variable with respect to 0 is zero. The origin of the coordinate system is located at the 
center of the accretion disk and the “z=0” plane is located to coincide with the accretion 
disk midplane. The hydrodynamic equations are presented in the Appendix. 

In this treatment we have implemented the energy equation self consistently. In partic- 
ular the adiabatic heating and cooling effects are implemented through the first and second 
term on the right hand side of equation (A6), and the radiative heating and cooling rates 
are implemented through the last term on the right hand side of equation (A6). 

As discussed in section 4, our treatment of the line radiation pressure assumes negligible 
the 4> component of the line radiation force, and therefore the specific angular momentum 
is conserved within our models. We note that strictly, the azimuthal component of the line 
radiation force is not required to be zero, and thus the angular momentum is not necessarily 
conserved, but this is likely a small correction. 

As we discovered in Paper 1 and Paper 2, care must be taken with the boundary 
conditions to avoid numerical instabilities. In this work we treat the boundary conditions 
in a form equivalent to that discussed in Paper 2. Our spatial grids are also equivalent to 
those explored and discussed in Paper 2. 

In the models presented here we solve the hydrodynamic equations over a range of radii 
from O.O2i?0 (2 R wd in our models) to 1.5 R e . From equation (6) it is found that the energy 
emission per area Q will increase slightly as radius decreases from 2 R wd down to * 1-36 R wd 
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and will decrease strongly as radius decreases from « 1.36R wd down to the white dwarf 
radii (R wd ). In Paper 1 we showed that the wind density tends to decrease with decreasing 
radiation flux. Thus the mass loss rate for radii less that 2 R wd will not alter significantly the 
total mass loss rate of the wind. Also, in Paper 2, we presented models with spatial grids 
with r down to R^j d and found essentially the same results with models with “r” down to 
2 Rwd, with model parameters otherwise identical. 


7. Results 

In Table 1 we present the wind mass loss rates from our previous models and for the 
models presented here. The second column indicates the dimension of the model, 2.5D are the 
full 3D hydrodynamic equations under azimuthal symmetry; i.e. the derivatives with respect 
to the azimuthal angle are assumed zero. 2.5D models allow a consistent treatment of the 3D 
equations (under azimuthal symmetry) while mathematically having only two independent 
spatial directions, keeping the calculations more manageable. The third column indicates 
the radial disk structure assumed. In all models a steady disk is assumed. An “I” indicates 
an isothermal disk (no radial dependence). An “S” indicates the radial dependence of disk 
radiation emission and temperature of a standard Shakura-Sunyaev disk. The third and 
fourth columns indicate the values of the k and a parameters used in each model. In our 
previous models (A-F) we had assumed constant ionization equilibrium throughout the wind, 
and therefore constant values for these two parameters. In our current work, as described 
above, we implement the spatial dependence (Model G) and the spatial and time dependence 
(Model H) of these two parameters and find values that vary throughout our spatial grids 
(and in time for Model H) within the intervals indicated in Table 1. 

In implementing local ionization equilibrium in our models, we first run models that 
include the spatial dependence of the radiation field above the disk but assume a fiducial 
number density of 10 8 cm -3 (Model G). From this we obtain the values of the “ k ” and “a” 
line radiation parameters for each point of the computational grid. Results from these runs 
are presented in Figure 5. The values of k we found varied from 0.03 to 0.5 and values of a 
varying from 0.6 to 0.8. The higher values of both k and a are obtained in the lower inner 
region of the wind (lower left hand side of Figure 5). The lower values of both k and a 
are obtained in the lower outer region of the wind (lower right hand side of Figure 5). At 
2 - R 0 , the k parameter was found to vary from 0.16 to 0.26 and the a parameter was found 
to vary from 0.66 to 0.70. In this model we find similar terminal velocities with respect to 
the models presented in Paper 2, but densities values increased by a factor of 2. We find 
steady wind flows coming from the disk. Wind-mass loss rates calculated in this work are 
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obtained by integrating the mass flux at the start of the wind, at the lower boundary of 
the spatial computational grid. The total wind mass loss rate we find for this model is 
2 x 10 -11 M q yr~T This increase in the wind mass loss rate is due to the increase of the k 
in the lower inner region where the radiation flux is stronger. 

In order to calculate the k{r, z, t) and a(r, z , t) functions self consistently throughout the 
hydrodynamic calculations, as we described in §5, the spatial distributions of the k and a 
are calculated in a similar manner as in the model presented in Figure 5, but now for values 
of number density varying from 10 6 cm- 3 to 10 10 cm' 3 at each computational spatial grid 
point (Model H). A fiducial electron temperature of 30, 000 K is taken for these calculations. 
During the hydrodynamic calculations, the actual values of k and a are obtained through 
geometric interpolation from the above values and the current value of density at each 
computational spatial grid point. Results from these models are presented in Figure 6 For 
the above range of densities at each spatial grid point, we find values of k varying from 0.03 
to 0.6 and values of a varying from 0.55 to 0.85. Again we find steady wind flows coming 
from the disk. In the lower inner region, due to the increase of density the k parameter 
drops by a factor of ~ 2, producing in turn a decrease of total wind mass rate. The wind 
mass loss rate for this model is 4 x 10~ 12 M e yr" 1 . The velocity fields are very similar to the 
previous models, finding terminal velocities in the order of ~ 3000 km s .We note that the 
mass loss rate is a approximately a factor of two below the values we found in our earlier 
models in Paper 2. This difference is due to the detailed calculations of the line radiation 
force parameters, which when the local density values at the inner region are considered, we 
find slightly lower values of the k line radiation force parameter. 

In Figures 7 and 8 we present contour plots of the line radiation force parameter k 
and a respectively, once the calculations arrive at a steady solution (Model H). These two 
figures illustrate the coupling between line force parameters, ionization equilibrium, and 
wind dynamics within an accretion disk wind. We note that, although the ionization struc- 
tures are coupled with the hydrodynamics, within typical CV parameters {M wd - O.6M 0 ; 
R wd = O.Oli?©; L disk = L q ), in the higher wind density areas, the line force parameters 
do not change significantly; thus justifying the assumption of our earlier models of con- 
stant ionization equilibrium (and thus constant line force parameters). We also note that 
in QSOs/AGN, the presence of strong x-ray emission will generate strong changes in t e 
ionization balance of the wind, and the assumption of constant ionization equilibrium in the 
wind will no longer hold (Murray et al. 1995; Hillier et al. 2002). 

Studies of observed ultraviolet line profiles from CVs have estimated CV wind mass loss 
rates to be in the order of 10~ 12 to 10- 11 M o yr- 1 and found CV wind terminal speeds in 
the order of ~ 3000kms _1 (e.g. Krautter et al. 1981; Cordova & Mason 1982; Greenstem 
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& Oke 1982; Mauche & Raymond 1987; Prinja & Rosen 1995; Friedjung, Selvelli, & 
Cassatella 1997) thus the models developed in this work finds values for wind mass loss 
rates and terminal speeds consistent with observations. 

In the models presented in this work we also implement radiative heating and cooling. 
In Figure 10 we present a velocity field graph superimposed with temperature contours for 
the local ionization equilibrium disk wind model. Comparing with the results from Paper 2 
(Figure 9), we find that radiative heating and cooling tends to maintain the wind temperature 
in the higher density regions near the disk, rather than cooling due to adiabatic expansion as 
we had found in our earlier models of Paper 2. Therefore the radiative heating and cooling 
processes are playing an important role in the temperature structure of the wind. We are 
finding higher overall temperatures than in Paper 2 where we neglected radiative heating 
and cooling processes and assumed only adiabatic heating cooling processes. 

As a consistency check, in Figure 11 we present the results of models with a computa- 
tional spatial grid that starts at r = R wd (rather than r = 2.07^ as in the other models 
discussed previously in this work). In Figure 12 we present results from models with a com- 
putational spatial grid starting at r = R wd and with an initial velocity at the base of the 
wind of lkms (rather than lOkms 1 as in the other models discussed previously in this 
work). In both cases we find similar overall results: a steady biconical wind velocity/density 
structure, with similar velocity and density values. We note that these results although 
similar, are not identical to those of Figure 6, and thus represent numerical effects since 
the physical parameters have remained identical. However the deviations from Figure 6 are 
considerably smaller than the changes of results we have found in earlier work when we 
physical parameters have been modified as our models have developed (Paper 2), thus we 
are confident in our overall results. In future work we shall study these numerical effects in 
greater detail; a precise knowledge of these model limitations may be important to achieve 
precise comparisons with observations. 

As we had found in our earlier models of Paper 1 and Paper 2, rotational forces are 
important in the study of winds from accretion disks. They cause the velocity streamlines 
to collide which results in an enhanced density region. In this region the wind speed is 
reduced and the wind density is increased. The increase in density caused by the collision of 
streamlines is important because it permits the appearance of blue-shifted absorption lines 
as observed in P Cygni profiles of low-inclination CVs. 
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8. Summary and Conclusions 

We have developed a 2.5-dimensional hydrodynamic line-driven accretion disk wind 
model. Our model solves a complete set of adiabatic hydrodynamic partial differential equa- 
tions, using the PPM numerical scheme and implementing the radial temperature and radi- 
ation emission distributions on the surface of an accretion disk. Our models calculate the 
line radiation force parameters throughout the wind at each computational grid point for 
each time step through the calculation of line opacities from atomic data, similar to the 
calculations presented by Abbott (1982) for OB stars. We treat the hydrodynamic equa- 
tion of energy self consistently including radiative heating and cooling as well as adiabatic 

expansion and compression. 

We find steady wind flows coming from the disk with similar density and velocity struc- 
tures as our earlier models of Paper 2 in which we had not yet included radiative heating 
and cooling and in which we had assumed constant ionization equilibrium (therefore constant 
values for the line radiation force parameters). The steady nature of our solution contrasts 
with the strong fluctuations in density and velocity found by Proga, Stone, k Drew ( 1998 ) 
for CV winds for physical parameters similar to the ones we used here. A possibility that, 
could account for the difference in results from our models and those of Proga, Stone, k 
Drew (1998) is the treatment of the lower boundary conditions. In future work we shall 

explore this possibility in detail. 

We note that as our numerical models evolve towards more realistic ones, it may re- 
sult that a more accurate treatment of the line radiation force (which would include line- 
overlapping, rather than assuming the Sobolev approximation for all lines) or a more accurate 
model of the accretion disk (which may result in an unsteady picture of the accretion disk 
rather than the standard Shakura-Sunyaev disk) may lead to physical instabilities. In our 
systematic approach to the accretion disk wind problem, we started with an analytic ID 
model, which although being our “crudest” one, did allow for an analytic solution which has 
proven to be invaluable as a test case for our numerical algorithms. From this initial model 
we have included, on a step by step basis (Paper 1, Pereyra 1997, Paper 2) elements such as 
a more accurate treatment of the line radiation force (in which we have always applied the 
Sobolev approximation), radial structure of the accretion disk, adiabatic heating and cooling 
processes, etc. This approach has allowed us to establish clearly a relationship between the 
elements we have introduced (both of numerical and physical nature) and our results, as 
our models have developed into more realistic ones. In this sense we are confident that if 
intrinsic physical instabilities should appear as we continue our work, we will be able to 
establish their physical origin. 

From our models we calculate wind mass-loss rates and terminal velocities. For typical 
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cataclysmic variable parameter values of disk luminosity L disk = L 0 , white dwarf mass 
M wd = 0 ,6M q , and white dwarf radii R wd = 0.01 Rq, we obtain a wind mass-loss rate of 
M wind - 4 x 10 12 Mq yr _1 , and a terminal velocity of ~ 3000kms _1 . We note that the 
mass loss rate is a approximately a factor of two below the values we found in our earlier 
models in Paper 2. This difference is due to the detailed calculations of the line radiation 
force parameters, which when the local density values at the inner region are considered, we 
find slightly lower values of the k line radiation force parameter. 

The results in this work suggest that the details of the ionization equilibrium in CV 
disk winds, within reasonable parameters, may not have considerable effect over the velocity 
fields of the wind nor in the wind mass distribution. Although these details have to be 
taken into account in order to derive an accurate theoretical value of the wind mass loss 
rate as illustrated by the difference by a factor of two between the models presented here 
and the earlier ones of Paper 2. The similar wind mass distribution and velocity fields found 
here with respect to those of Paper 2, would indicate a similar inclination angle dependence 
of the C IV 1550 A line. But the decrease in total wind mass loss rate would generate 
weaker line features. In QSOs, where the ionization balance of the gas material may change 
significantly, the coupling between the ion populations and the radiative hydrodynamics will 
have a important effects in the overall wind dynamics; thus a model capable of representing 
this coupling becomes necessary. 

A shortcoming of our models is that, at each spatial point for each time step we assume 
direction-independent values for the line radiation parameters and we also assume position 
independent values for the radiative heating and cooling parameters. In future work we plan 
to include these effects. We also plan to extend our disk wind models to low mass X-ray 

binaries and QSOs/AGN, where the local ionization equilibrium will play a more crucial role 
in the overall dynamics. 
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A. Hydrodynamic Equations 

The hydrodynamic equations implemented in the models presented are, the equation of 
state: 


p = ( 7 - 1 )pe 


(Al) 
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the mass conservation equation: 


dt r dr 

the momentum conservation equations: 
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and the energy conservation equation: 
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+ 


aS z (r,z) 


k(r, z, t) ^max 


1 

.PVVth 


dvz 

dz 



a(r,z,t) 

V z 


+ n 2 (r(T) — A(X’)) 

where P is the pressure, 7 is the ratio of specific heats, p is the density, e is the internal energy 
per mass, v T , v$, and v z are the corresponding velocity components in cylindrical coordinates, 
G is the gravitational constant, M wd is the mass of the white dwarf, a is the Thomson cross 
section per unit mass, F r and F z are the corresponding radiation flux components (see 
equations [42] and [43]), c is the speed of light, v t h is the thermal velocity, and k and a are 
the line force multiplier parameters calculated through equation (30), S r and S z are defined 
through equations (38) and (37) respectively, E = v 2 r /2 + n 2 /2 + v 2 z /2 + e is the total energy 

per mass, n is the electron number density, T is the radiative heating rate, and A is the 
radiative cooling rate. 
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Table 1. CV Accretion Disk Wind Mass-Loss Rates for Different Models 


Model 

Dim 

Disk 

k 

o 

A 

1 

I 

1/30 

0.7 

B 

2.5 

1 

1/30 

0.7 

C 

2.5 

I 

1/3 

0.7 

D 

2.5 

I 

1/3 

0.7 

E 

2.5 

I 

1/3 

0.7 

F 

2.5 

S 

1/3 

0.7 

G 

2.5 

S 

0.03-0.5 

0.6-0. 8 

H 

2.5 

s 

0.03-0.6 

0.55-0.85 


Comments 


Isothermal Wind 
Model A extended to 2.5D; 
radiation pressure in “r” neglected 
Model B with an increased value of k 
Model C with a maximum cutoff 
value for the line pressure implemented 
Model D with an adiabatic wind 
(rather than an isothermal wind) 

Model E with a 

standard Shakura-Sunyaev disk 
(rather than isothermal disk) 
radiation pressure in “r” included 
Model F with spatial 
dependence of k and a derived 
from the spatial distribution 
of radiation field but assuming 
a fiducial number density of 10 s cm 3 
Model F with spatial and time 
dependence of k and a derived from 
the spatial distribution of radiation and the 
spatial and time distribution of mass density 


M(Mq yr~ l ) Reference 


2 x 10 -14 
2 x 10 -14 

Paper 1 
Paper 1 

6 x 10~ 14 
2 x 10~ 13 

Paper 2 
Paper 2 

2 X 10~ 13 

Paper 2 

8 x 10~ 12 

Paper 2 

2 x 10~ u 

current work 

4 x 10~ 12 

current work 
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log v (Hz) 


Logarithm of radiation flux vs logarithm of frequency in units of Hz for the accretion 
disk used in this work. This is the spectrum distribution used for the radiative heating and 
cooling calculations. The physical parameters used here are L disk = L e and R wd = 0.01i? o . 
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Pig 4. Vector field graph of wind velocity superimposed with density contours for the 

model assuming constant ionization equilibrium (Model F [Paper 2]). The line radiation 
force parameters applied are k = 0.3, a = 0.7 . The primary star is at the origin of the 
graph, and the disk is over the horizontal axis. The contour levels vary uniformly from a 
value of 3.6 x 10 -16 gcmT 3 down to a value of 0.1 x 10 -16 gcm -3 . The physical parameters 
here used are M n ,d — O.6M0, R w d — 0.017?©, and L^isk — T©. 
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Fig. 5. Vector field graph of wind velocity superimposed with density contours for the 
model assuming an ionization equilibrium varying in space but constant in time. A fiducial 
number density of 10 8 cm 3 is used for the calculations of the line radiation force parameters 
(Model G). Parameters are similar to those of Figure 4. 
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Fig 5 Vector field graph of wind velocity superimposed with density contours for the 

model applying an ionization equilibrium varying in space and time. The current density is 
applied at each spatial grid point at each time step in the calculation of the line radiation 
force parameters (Model H). Parameters are similar to those of Figure 4. 
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pig 8.— Contour plot of the line radiation force parameter a in the Model H once it arrives 
at a steady solution. The physical parameters used are those of Figure 4. 
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Fig. 9. Vector field graph of wind velocity superimposed with temperature contours for the 
model assuming constant ionization equilibrium and without radiative heating and cooling 
(Model F [Paper 2]). The contour levels vary uniformly from a value of 22, 000 K down to a 
value of 1, 000 K. Parameters are similar to those of Figure 4. 




r (R e ) 

Fig. 10 — Vector field graph of wind velocity superimposed with temperature contours for 
the model applying local ionization equilibrium and including radiative heating and cooling 
(Model H). Parameters are similar to those of Figure 9. 
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Fig. 11. Vector field graph of wind velocity superimposed with density contours for the 
model applying local ionization equilibrium varying in space and time. The computational 
grid starts at r = rather than r = 2R wd as in the results presented in Figure 6. 

Parameters are similar to those of Figure 4. 
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Pig 12. — Vector field graph of wind velocity superimposed with density contours for model 
disk wind applying local ionization equilibrium varying in space and time. The computational 
grid starts at r = R wd , rather than r = 2 R wd as in the results presented in Figure 6. The 
initial velocity set at the base of the wind is 1 km s \ rather than lOkms as in the results 
presented in Figure 6. Parameters are similar to those of Figure 4. 




